Doping the Kane-Mele-Hubbard model: A Slave-Boson Approach 
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We study the Kane-Mele-Hubbard model both at half-filling and away from half-filling using a 
slave-boson mean-field approach at zero temperature. We obtain a phase diagram at half-filling 
and discuss its connection to recent results from quantum Monte Carlo, cellular dynamical mean 
field, slave-rotor, and Z2 mean-field studies. In particular, we find a small window in parameter 
space where a spin liquid phase with gapped spin and charge excitations reside. Upon doping, we 
show the spin liquid state becomes a superconducting state by explicitly calculating the singlet 
pairing order parameters. Interestingly, we find an "optimal" doping for such superconductivity. 
Our work reveals some of the phenomenology associated with doping an interacting system with 
strong spin-orbit coupling and intermediate strength electron-electron interactions. 
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I. INTRODUCTION 

Recent years have seen growing interest in topological 
band insulators (TBI). 1-3 While in the fractional quan- 
tum Hall effect the Coulomb interaction is necessary for 
the topological order, the concept of TBIs can be read- 
ily appreciated in the framework of noninteracting Bloch 
band theory where spin-orbit coupling is responsible for 
a possibly nontrivial Z2 topological order. 1_3 A TBI has 
a gap in the bulk excitation spectrum and time-reversal 
symmetry protected gapless modes on the boundary. Be- 
cause in nature all materials possess Coulomb interac- 
tions, understanding the role of interactions is one of the 
fundamental challenges in the held of topological insula- 
tors. 

One can ask if it is possible for interactions to induce 
TBIs. The answer is affirmative. Indeed, there have 
been a number of works that address this question with 
different models that contain no intrinsic spin-orbit cou- 
pling. For example, Raghu et al. showed it is possi- 
ble to have an interaction-driven TBI with spontaneously 
broken SU(2) symmetry (with spontaneously generated 
spin-orbit coupling) from an extended Hubbard model on 
the honeycomb lattice. This idea has been successfully 
applied to the kagome lattice 5 and the decorated honey- 
comb lattice 5 in 2D among others, 6 ~ 12 and the diamond 
lattice in 3D. 7 The key is to have the correct amount 
of "generalized" spin-orbit coupling that originates from 
the Hartree-Fock mean- held decoupling of the interaction 
terms on nearby sites. 

Another equally important question is the fate of TBIs 
with intrinsic spin-orbit coupling upon the inclusion of 
Coulomb interaction. On one hand, by the argument 
of adiabatic continuity, it is argued that a TBI should 
be stable to weak interactions as long as the bulk gap 
is not closed. 13,14 However, when interactions grow too 
strong, one has a good reason to believe that spin-charge 
separation develops and Mott physics will appear. 15 In 
this regime, one expects that a slave particle approach 



which starts with an explicit decomposition of the elec- 
tron into charge and spin degrees of freedom would qual- 
itatively capture the physics of the interactions. Indeed, 
back in 2008 Young et al. 16 employed a slave-rotor mean- 
field approach to study a double layer honeycomb lattice 
where a fractionalized quantum spin Hall (FQSH) effect 
could be found. A FQSH state differs from a quantum 
spin Hall state in that neutral spinons instead of physical 
electrons carry a nontrivial Zi topology. As a result, a 
gapless spinon excitation is guaranteed to appear along 
the edge. Applying similar methods, others 17,18 used the 
same approach to study the Kane-Mele-Hubbard model 
[our Eq.(l)] on the single-layer honeycomb lattice and 
concluded that this phase could be stabilized if the two 
dimensional U(l) gauge field is screened by an additional 
metallic layer so that the gauge fluctuations are sup- 
pressed. 

In three dimensional systems, Pesin and Balents 15 
studied heavy transition-metal oxides on the pyrochlorc 
lattice and proposed a three dimensional counterpart 
of the FQSH, termed as a "topological Mott insula- 
tor" (TMI). A TMI is one example of a U(l) spin liq- 
uid (SL) in three dimension and is believed to be more 
stable to gauge fluctuations than its two dimensional 
counterpart. 19 Later, Kargarian et al. 20 extended Pesin 
and Balent's results and investigated the interplay be- 
tween interactions and distortion in the same system. 
Based on these works, it may appear that the concept 
of the FQSH in two dimensions and the TMI in three di- 
mensions depends crucially on the slave-rotor approach, 
which by its construction transfers the topology of physi- 
cal electrons to neutral spinons and makes access to frac- 
tionalized states possible. 

The extent to which a slave-rotor mean-field approach 
is reliable can be checked with more controlled numeri- 
cal simulations. Recent quantum Monte Carlo and cel- 
lular dynamical mean field studies have shed light on 
the weak and intermediate interaction regimes in two 
dimensions. 21 ~ 25 In a pioneering quantum Monte Carlo 



study, Meng et al. 21 investigated the Hubbard model on 
the honeycomb lattice at half-filling and discovered the 
existence of a gapped spin liquid in a small window in 
the intermediate interaction regime (3.5i < U < 4.3i). 
Later, spin-orbit coupling was included and the spin liq- 
uid phase was found to be stable for small spin-orbit 
coupling 22 and for finite temperatures. 25 At half-filling, 
the above quantum Monte Carlo studies are free of the 
sign problem and considered to be accurate. Of par- 
ticular interest is the nature of the spin liquid, which 
has been addressed in a number of works. 26 ~ 30 Very re- 
cent work has indicated that beyond a critical interaction 
strength and spin-orbit coupling strength (larger than 
that explored in quantum Monte Carlo) yet another novel 
phase may appear with fractionized excitations and a 
non-trivial ground-state degeneracy 31 and attention has 
been drawn to transition metal oxide interfaces. 10 ~ 12,32 



II. THE SLAVE-BOSON APPROACH 

We start with the Kane-Mele-Hubbard model on the 
honeycomb lattice, 



H 



* 2-^1 C i..a C i- (T 



■U 






n^riii- 



-iX 



so 



^ cri^-ct c. 



«> 



where t, U, and Aso are the nearest neighbor hop- 
ping energy, the strength of the on-site repulsion, and 
the second-neighbor spin-orbit coupling strength, respec- 
tively. Here Cj CT (c ia ) annihilates (creates) an electron 
with spin a on site i and i>ij — ±1 depending on if the 
electron makes as "right" or "left" turn when going from 

i to J. 38 ' 39 

The general U(l) slave-boson approach decomposes an 
electron operator into a bosonic operator that carries the 
charge degree of freedom and a fermionic spinon operator 
that carries the spin degree of freedom: 28 ' 36,40-42 
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In this paper we aim to better understand the inter- 
mediate interaction regime where a gapped SL phase 
appears. We are particularly interested in the fate of 
the SL 33 upon doping. This is a regime where quantum 
Monte Carlo simulations suffer from the sign problem 34 
and the slave-rotor mean-field approach may encounter 
severe limitations 35 leaving few tools available for its 
study. We will follow Ref. [28] and use a generalized 
U(l) slave-boson mean-field approach to study the cases 
of half-filling and doping. Such an approach has been 
widely used in doped t-J models in the context of high 
temperature superconductivity. 36 We stress that we do 
not expect the slave-boson mean-field approach to repre- 
sent a good solution to the Kane-Mele-Hubbard model in 
all regimes. Instead, we argue that it gives a reason- 
ably good description of the gapped SL at intermediate 
regime (based on a quantitative comparison with QMC 
and CMDFT) and its transition to a superconducting 
state upon doping. For a general review of Hubbard 
model, we refer interested readers to Ref. [37]. 



This paper is organized as follows. In Sec. II we intro- 
duce the slave-boson representation for the Kane-Mele- 
Hubbard model. In Sec. Ill we describe our slave-boson 
mean-field results for the cases of half-filling and dop- 
ing. Finally, in Sec. IV we give the main conclusions 
of this work. In App. A we provide some lengthy self- 
consistency formulas used to obtain our results. 



where i is the site index, and hi and di are the bosonic 
holon operator and the bosonic doublon operator, respec- 
tively. Such a decomposition makes the idea of spin- 
charge separation explicit and one expects that it will 
describe the physics of intermediate (and possibly strong) 
interactions reasonably well. 

There are four states, |0),| f)>| \) 7 an( l I t4-)s at eacn 
site. Each state can be thought to have some new parti- 
cle operator acting on some vacuum state: |0) = h'\vac), 
I t) = fl\vac), I |) = fl\vac) and | U) = <$\vac). Physi- 
cally, one can think of h\hi as the number of empty oc- 
cupancies at site i, fj^fia- the single occupancy with spin 
a, and d\di the double occupancy. One can show that 
Eq. (2) guarantees that the matrix elements of physical 
states are correct. The completeness of the basis implies 
the constraint 
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which also preserves the anticommunication relations of 
Ci jCr and c ia . This constrain can be enforced with a La- 
grange multiplier A^ in the Hamiltonian. 

There is also another constraint related to the filling 
fraction of electrons: c ia Ci a — f i(J fia + d\di where some 
extra terms which have zero matrix elements in the phys- 
ical states have been thrown away. 28,36 ' 40-42 Therefore, 
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where x is the electron doping. We can incorporate the 
constraint (4) by another Lagrange multiplier \ii in the 
Hamiltonian. 

With the slave-boson representation described above, 
the Kanc-Mclc- Hubbard model can be written as 
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where the following order parameters are defined for 
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We proceed with a mean-field approximation in which 



the spinon part and the boson part decouple from each 
other. We will restrict ourselves to a search for phases 
preserving translational symmetry. The simplest phase is 
the one that does not break any symmetry of the Hamil- 
tonian. In this case each type of order parameter does 
not depend on the site indices. For example, (xL) = Xb 



for any (ij). 

where 



We then have H MF = H f + H b + H c . 



H f = E \_- t Xbg{^)ft AlJ hBa - iA h g(k)o-/^ Aff /l kB _ ff + \soXb<r(fkA*faAa - fl B afaBa)gi(k) 

Via 

a E fL*fi™ + Aso a^ E f/iLt/Ik^ - fUf-±Bi) + h - c -> 



(6) 



kcra 



Hb = ^2\{U - \- fi- A S oX/52(k)) dl a d ka + (/i - A + A S oX/52(k)) hl a h ka + A S oAy5 2 (k)d kiQ /i_ kiQ 



E Vx}{h\ A hkB - d kj4 rfks)g(k) + tA/(d k Aft-kB + ZikAd-ksM-k) 



/i.e. 



(7) 



r 



and the constant energy term H const = 6Nt(xtXf + 
A fc A/) + 2N{X + xij,) ~ 12N( x ' b x'f + A 'b A 'f) wh ere N 
is the number of unit cells. In Eq.(6) and (7) we have 
used a = A, B to denote the two sublattices of the hon- 
eycomb lattice. We have defined g(\s) = 1 + cxp(— ik%) + 
exp(zfci — ifc 2 ), .9i (k) = 2[sinfc 2 — sinfci+sin(fci — fe)] and 
g 2 (k) = 2 [cosfci + cosfc 2 + cos(/c 1 — fe 2 )] vvith fcj = k^. 
We stress that both spinon and bosonic Hamiltonians 
have generic hopping terms and pairing terms and they 
are not identical to the pairing of physical electrons, as 
we will explain later in our paper. 

After a mean-field approximation, it is straightforward 
to solve the spinon Hamiltonian (6) and the bosonic 
Hamiltonian (7) to obtain the ground state energy at 
zero temperature. Self-consistency equations are ob- 
tained via the first derivative of the ground state energy 
with respect to various order parameters. (See the ap- 
pendix for details.) However, one has to consider possible 
Bose- Einstein condensations when dealing with a bosonic 
Hamiltonian. Since we have assumed that the translation 
symmetry remains unbroken, we expect a Bose-Einstein 
condensation can only take place at k = 0. Bearing this 
in mind, we can explicitly separate the k = term from 
the k / terms. We define d a = l/vN(dk=o,a) an d 



h a = 1/v N(hk—o, a ) j both of which can acquire finite 
values in a condensed phase of the bosons. Thus, d? a and 
h 2 a are the fraction of doublons and holons on sublattice 
a = A,B. 

We will also consider phases that break certain symme- 
tries (lattice rotational symmetry, for example) at large 
U. This allows us to make connections to an antifer- 
romagnetic state, which is difficult to capture within 
the slave-boson mean-field approach to the Kane-Mele- 
Hubbard model. 28 In the next section we turn to a de- 
tailed description of the results of our mean-field study. 
We find many features reminiscent of previous studies at 
half-filling, but we also obtain new results for the doped 



III. MEAN-FIELD RESULTS 



In this section we discuss our mean-field results for the 
cases of half-filling and doping away from half-filling. We 
begin with the half-filled case. 
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FIG. 1. Phase diagram of Kane-Mele-Hubbard model at 
half filling within the slave-boson mean-field approach. SC 
stands for the superconducting phase, SL is the spin liq- 
uid phase and DM is the dimerized phase. The DM phase 
may be considered as the closest relative to the spin den- 
sity wave (SDW)/antiferromagnetic state obtained in pre- 
vious studies. 17 ' 21 ~ 25 Importantly, the slave-boson treatment 
also produces a gapped spin liquid at intermediate coupling 
(which was found in earlier numerical studies 21-25 ). We note, 
however, that the slave-boson treatment does not smoothly 
connect to the non-interacting limit since it predicts a SC 
phase rather than a TBI. This is a shortcoming of the slave- 
boson mean-field method which is designed to capture the 
physics of intermediate U values where the SL phase appears. 



A. Half-filling case 

For the case of half-filling our results are summarized 
in Fig. 1. There are many similarities with results ob- 
tained in the literature via different techniques. 17 > 21 ~ 25 
Most importantly, we find a gapped spin-liquid phase at 
intermediate coupling that extends to finite spin-orbit 
coupling. In order of increasing interactions, the phases 
we find are: 

(1) Superconducting states (SC)-When the interaction 
strength U is small, there is a finite probability of dou- 
ble occupancy and empty occupancy at each site. There- 
fore, we expect Bose-Einstein condensation of holons and 
doublons could take place for small U. Indeed, we find 
a critical interaction strength U c (Xso) an d U c ~ 1.5* at 
Aso = above which the SC phase does not survive, 
which is about half of the value that has been reported 
in the quantum Monte Carlo simulation. 21 (Although in 
that case it is a semi-metal that persists up to a critical 
interaction strength.) The SC phase persists even for a 
negative interaction, as one might expect. One can view 
the U > SC region as an "extension" from the U < 
region to "small" repulsive interactions. Recent argu- 
ments have shown that SC can, surprisingly, be expected 
even for (small) repulsive interactions. 43 However, as we 



emphasized earlier, the slave-boson mean-field treatment 
does not properly capture the small U > physics prop- 
erly in the model at half-filling. We do expect SC states 
to be likely for small U upon doping, and for that reason 
also discuss the technical details of the half-filled case 
here which will only be slightly modified upon doping. 

From mean-field self-consistency equations, we find 
that this superconducting state can be described by four 
finite condensates Ha, hs, d,A and ds and finite A/, At>, 
A'f and A' b (SC I). The four condensates are related via 
Ha = —hs and <i^ = — ds (or other equivalent config- 
urations). All other order parameters (i.e. the \) are 
zero. The physical picture for small interaction is then 
as follows: the spinons are paired at nearest and sec- 
ond nearest sites and cannot hop freely on the lattice; 
bosons (doublons and holons) condense independently at 
k = in momentum space. The ground state has gap- 
less charge excitations and gapped spinon excitations. In 
terms of the physical electrons pairing, one can show that 
generally, 
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To further discuss the properties of this SC phase, let's 
consider singlet pairing between the same sublattices in 
the absence of spin-orbit coupling and we find 



C 0at C r«V 



-Ae 



— ik-j 



k 



\/A 2 + t 2 | ff | 2 A 



2' 



(9) 



where we have used fact that the Bose-Einstein conden- 
sation takes places at weak interactions so that we can 
replace bosonic operators with their averages. Therefore, 
one has finite on-site and NNN singlet pairings between 
same sublattices and also for neighbors arbitrarily far 
away. On the other hand, we find singlet pairings be- 
tween different sublattices vanish. It is also possible to 
obtain another SC solution with finite %s and conden- 
sates but zero As (SC II). 28 The spinon sector is the 
effective noninteracting Kane-Mele model with physical 
electron operators replaced by neutral spinon operators. 
Clearly this spinon Hamiltonian possesses non-trivial Z 2 
topology and has time-reversal symmetry protected gap- 
less edge states. However, singlet pairings for electrons 
between same sublattices are finite, 
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where E s (k) = —A + stx&|<7(k)| and O is the Heaviside 
step function. Similar to SC I, SC II has zero singlet 
pairings between different sublattices. That's the reason 
we identify it as a SC state. However, we find SC II is not 
energetically favorable. In Fig. 2, we explicitly show the 
difference of two mean-field solutions. Note: our mean- 
field solutions at half filling only admit the above two 
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FIG. 2. (Color online) The ground state energies for two 
slave-boson mean- field solutions. 



solutions and there exists no phase with ^^0 and A ^ 
at half filling. 

As pointed out in Ref. [28], in the weak interacting 
limit the Bose gas of doublons and holons is dense (large 
amplitude of condensates) and one should expect the ex- 
istence of strong interactions between them. Therefore, 
the slave-boson mean-field approach is not reliable for 
weak interactions. Indeed, the ground state for weak 
interactions in the absence (presence) of spin-orbit cou- 
pling is a Fermi liquid (TBI) . This is confirmed in a recent 
quantum Monte Carlo study. 22 

Another popular approach to handle interactions, the 
slave-rotor mean-field approach, is believed to be able to 
reasonably capture the qualitative features of physics at 
small interactions. 35 ' 44 It has been applied to the Hub- 
bard model on the honeycomb lattice and predicts a 
nodal spin liquid phase for 1.68t < U < 1.74£. 45 Later, it 
was applied to Kane-Mele-Hubbard model on the same 
lattice and successfully predicted a TBI phase for weak 
interactions, though the gauge field has to be screened 
out to stabilize it. 17 The mathematical structure of slave- 
rotor approach allows a direct transfer of topology from 
physical electron bands to neutral spinons; this is the key 
to predicting a TBI at weak interactions and a TMI at 
intermediate to strong interactions. However, the slave- 
rotor method suffers from severe limitations for finite 
doping. 35 The slave-boson mean-field approach, on the 
other hand, allows in principle nontrivial band topol- 
ogy embedded in its spinon sector. Unfortunately, in our 
case we obtain only finite paring terms. As a result, the 
slave-boson mean-field approximation falsely predicts a 
SC for half-filling and weak interactions. We also want to 
mention the Kotliar-Ruckenstein slave-boson mean-field 
approach describes the weak interacting limit well, 46,47 
though it might be difficult to address the intermedi- 
ate coupling regime and obtain a gapped spin liquid. It 
would be interesting to study its predictions and this will 



FIG. 3. Double occupancy as a function of U for A = Q.Q2t 
and half-filling at zero temperature from slave-boson mean- 
field approach. Note the double occupancy at the weak inter- 
acting limit is not correct since it is larger than 1/4, but it 
gets better as the interaction grows. 



be left as a future work. 

(2) Spin liquid states (SL)-As the interaction grows, 
the Bose gas becomes less dense, and one expects that 
the slave-boson approach is better able to describe the 
intermediate interaction regime. We find a spin liquid 
phase appears between 1.5i < U < 1.9t for Aso = 0. 
In the absence of spin-orbit coupling, this phase is char- 
acterized by finite A/ and A&. Both the spinon sector 
and the chargeon sector are gapped and no Bose-Einstein 
condensation takes place. We find the singlet pairings be- 
tween any two sites vanish. The expectation value of the 
spin at each site is also zero, and the spin-spin correlation 
decays exponentially due to a finite spinon gap. There- 
fore, we obtain a spin liquid phase in a small interaction 
window. Furthermore, we find it can survive over a small 
range of spin-orbit coupling. This feature is quite similar 
(even numerically) to the quantum Monte Carlo result, 
though the specific phase boundary differs. 22 

To substantiate our assertion that the slave-boson 
mean-field approach gets better when the interaction 
grows, we follow Ref. [25] and plot the double occupancy 
D occ = (ni-fUii) for Aso = 0.02t and half-filling at zero 
temperature in Fig. 3. As one can see, in the weak in- 
teracting regime, D occ is larger than 1/4 (which is the 
exact values for U — 0) and this is another evidence that 
slave-boson mean-field approach does not work well in 
the weak interacting regime. However, as the interac- 
tion grows, for example, at U = 1.9i, our D occ — 0.23 
at zero temperature and this can be compared with Ref. 
[25] 's D occ « 0.21 at T = 0.025i. Since a finite temper- 
ature tends to reduce the double occupancy, we expect 
that our result will be very close to that of Ref. [25] if a 
zero-temperature cellular dynamical mean field study is 
performed. 



To study the SL phase in more detail, we calculate 
the single particle retarded Green's function G r aa (k, r) = 
— i6(t){{ck aa (t),c^ aa }) in the absence of spin-orbit cou- 
pling for the SL phase and the result is 
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where E(q, k) = £^(q) + _E fc (q-k) contains a ferminonic 
excitation _E^ (k) = ^/A 2 + i 2 |g(k)| 2 A^ and a bosonic ex- 
citation E b (k) = J(U/2- \y 
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|(1 + E f X , k s ) ■ As we are considering half-filling, the re- 
tarded Green's function exhibits particle-hole symmetry. 
To make our calculations more solid, we first check if 
the sum rule of spectrum function p(k,u;) = -Im[G^ CT ] 
is satisfied. Since it is based on the anticommunication 
relations between Ck and c k and it has been taken into 
account by Eq. 3 on the average, the sum rule of our 
slave-boson mean-field approach is implicitly fulfilled by 
the mean-field equations. The local density of states 
p(cu) = X)k' 9 (' t ' a; ) 1S snown m Fig- 4 for U = 1.8t and 
-^so = 0. The single particle gap is found to be 0.57t. 
This is the gap at the Dirac point. Instead of calculating 
it numerically in Ref. 21, we can determine it analytically 
in our case. The poles of the retarded Green's function 
are at w = ±£'(q, k) and the positive pole reaches its 
minimum at Dirac point, q = k = K, therefore it is 
clear that the single particle gap at the Dirac point is 

A sp = E(K,K) = |A| + J(U/2 - A) 2 - 9i 2 A 2 = 0.57* 



for U — 1.8t. Since our phase boundary for SL differs 
from Ref. [21], we cannot compare the single particle gap 
directly for the same U. However, our result for a point 
sitting about in the middle of SL phase (0.57i) is compa- 
rable to a typical single particle gap from Ref. [21] (about 
(Hi). 

A SL is also found in the slave-rotor mean-field ap- 
proach, though its properties are quite different from 
the one obtained here. 17,18 For example, only hopping 
terms of spinons are present in the SL within the slave- 
rotor approach and it has a U(l) gauge symmetry. In 
2D, U(l) gauge fluctuations are important 19 and it has 
been argued that other gapless layers (spinons) are re- 
quired to screen the gauge field and suppress the gauge 
fluctuations. 16 Our spin/charge gapped SL, however, 
does not require an additional layer to stabilize the mean- 
field result. The key difference with the slave-boson ap- 
proach is that the effective fermionic Hamiltonian for 
the SL consists of pairing terms instead of hopping of 
spinons. As a result, the presence of NNN pairings al- 
low the staggered U(l) gauge symmetry to break down 
to a Z 2 gauge symmetry by the Anderson-Higgs mecha- 
nism and gap out the U(l) fluctuations. 28 Therefore, our 
mean-field result has at least a chance of being realistic, 
and quantum Monte Carlo calculations 22 and dynamical 
mean-field theory calculations 25 support this result in a 




FIG. 4. The local density of states for the spin liquid phase 
at half-filling. We have taken U = 1.8t and Aso = 0. 



similar parameter regime. 

(3) Dimerized phase (DM)-The spin liquid phase is 
unstable to dimcrization in the large U limit. The 
dimerized phase is close in many respects to an an- 
tifcrromagnetic phase, which is expected to occur at 
large interactions on a bi-partite lattice like the hon- 
eycomb lattice. 17 > 21 ~ 25 With the present form of Kane- 
Mele-Hubbard model (in the absence of a spin-exchange 
term), it is difficult to include antiferromagnetic order 
in our mean-field approach. 28 We will instead turn to a 
dimerized phase which has anisotropy in some direction 
(i.e. breaks lattice rotational symmetry) yet keeps the 
translational symmetry intact. Similar ideas have been 
applied in the slave-rotor approach. 16 To perform our cal- 
culations, we will assume rotational symmetry is sponta- 
neously broken. We consider an ansatz of three different 
mean-field A/i, A/2, and A/3 for NN pairings, and A'^, 
A/2, and A'/ 3 for NNN pairings. We do the same in 
the chargeon sector. In the parameter space we consider, 
the mean-field solutions are those that satisfy A/i 7^ 0, 
A/ 2 = A/3 = and A bl ^ 0, A h2 = A h3 = while 
the NNN pairings for spinons and bosons vanish. This is 
an extreme example of dimerization and it corresponds 
to an atomic-like insulator which consists of noninteract- 
ing pairs of NN sites. This can be taken as a "closest 
cousin" to the antiferromagnetic state expected at large 
U for half-filling. 

To make further connections to the numerical studies, 
we follow Ref [21] and plot the derivative of the kinetic 
energy per unit cell dEkin/dU as a function of U/t in 
Fig. 5. After a comparison to the QMC, we find: (i) our 
kinetic energy is higher than the one in QMC and we 
expect our ground state energy is also higher, though we 
are not aware of reported ground state energy in QMC; 
(ii) Our kinetic energy profile resembles the one in QMC, 
though we have a jump around U = 1.9i from the spin 
liquid phase to the dimer phase while it has a continuous 
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FIG. 5. The derivative of kinetic energy per unit cell 
dEkin/dU for Aso = 0. Insert: the kinetic energy and the 
ground state energy. The kinetic energy has a jump around 
U — 1.9* and the derivative of it shows a sharp peak at the 
same location, which indicates a first order transition around 
U = 1.9*. 



behavior in QMC; (iii) Our calculations show that we 
have a second order phase transition at the first critical 
point U c i = 1.5* followed by a first order phase transition 
at U C 2 = 1.94 while there appears to be a continuous Mott 
transition around U = 3.5i. 26 



B. Doping Cases 

Up to this point, we have focused on the case of half- 
filling and our mean-field results could be directly com- 
pared with other numerical approaches. 17 > 21 ~ 25 We now 
break new ground by considering the case of doping away 
from half-filling where other methods may encounter se- 
rious shortcomings. 

The doped Hubbard model in the strongly interact- 
ing limit and its descendant t-J model are believed to 
capture the physics of high temperature superconduc- 
tivity upon doping. 36 ' 48 In most slave-boson treatments, 
one assumes strong interactions and simplifies the calcu- 
lations by removing double occupancy from the Hilbert 
space. However, since we are mostly interested in the in- 
termediate regime where U and * are comparable, and a 
possible spin liquid phase resides, we will start with the 
Kane-Mele-Hubbard model without assuming a strong 
interaction. Wc therefore retain the entire Hilbert space. 

Using the mean-field self-consistency equations (see 
Appendix for details), one finds that the SL at half-filling 
is unstable to infinitesimal doping and a Bose-Einstein 
condensations of chargeons takes place for any doping. 
This can be seen from Eq. (A6) where the doping is di- 
rectly related to the condensates. The number of dou- 
blons at each site is not equal to the number of holons, 



FIG. 6. The NN order parameters as functions of doping. 
We have set U — 1.8* and Aso = 0.05*, which is a SL at 
half-filling. 



and at least one of them has to be finite. This indicates 
the onset of Bose-Einstein condensation for any doping. 
Our mean-field solutions show that the \ s also acquire 
finite values, i.e. spinons and chargeons can both hop 
and form pairs on the lattice. 

In Fig. 6 and Fig. 7, we show various NN and NNN 
order parameters. As one can see, Xb and x'b have linear 
relations with respect to doping, which readily follows 
from Eq. (A6), Eq. (A7) and Eq. (A8). Xf and x'f have 
similar behaviors and are odd functions of doping while 
Af,, Af,A' b and A* are even functions of doping. Inter- 
estingly, the value of A' b and A', are numerically very 
close to zero at half filling. The four condensates are re- 
lated via hj\ = —hs and d^ = -de (or other equivalent 
configurations) . 

In Fig. 8, we plot the physical onsite, NN and NNN 
singlet pairings as a function of doping for parameters 
U = 1.8* and Aso = 0.05*, whose ground state is a 
SL without doping. As the doping is increased, singlet 
pairings between same sublattices and different sublat- 
tices acquire finite values and signal the occurrence of 
a SC phase. The singlet pairings are not monotonic 
functions of the doping and there exists an "optimal" 
doping (around ±0.8 electron/site) where the magnitude 
of on-site and NN parings are maximized. This bears 
some similarity to the famous SC "dome" in the phase 
diagram of high temperature superconductors, 36 though 
electron doping and hole doping are "equivalent" in our 
case. We also remark that the dimcrized phase will be- 
come a SC state via doping. Therefore, upon doping 
the SC phase takes over the entire phase diagram within 
the slave-boson mean-field treatment. However, as we 
mentioned earlier, the SC phase obtained is not one that 
possess topological order of any obvious type. 

In Fig. 9, we plot the ground state energy as a func- 
tion of doping. We have set U = 1.8* and Aso = 0.05*, 
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FIG. 7. The NNN order parameters as functions of doping. 
We have set U — 1.8£ and Aso = 0.05i, which is a SL at 
half-filling. 



FIG. 9. The ground state energy per unit cell E g as a function 
of doping. We have set U = 1.8t and Aso = 0.054, which is a 
SL at half-filling. 
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FIG. 8. (Color online) Singlet parings at U = 1.8f and Aso = 
0.05£ corresponding to the SL at half-filling in Fig.l. Shown 
as a function of doping (additional electrons/site) is: (a) The 
on-site pairing, (b) the nearest neighbor paring, and (c) the 
next nearest neighbor pairing. The black solid line is for the 
real part of the pairing, which is identical for both sublattices; 
the blue and red solid (dash) lines are for the real (imaginary) 
part of paring that is different for A and B sublattices. 



which is a SL at half-filling. For the ground state en- 
ergy, one can understand it as follows. At x = — 1 where 
electrons are completely depleted the energy is zero, and 
when one starts to add more electrons to the system, the 
ground state energy decreases since kinetic energy dom- 
inates over the potential energy and lowers the ground 
state energy. As more electrons are added, the poten- 
tial energy starts to dominate and cause the increase of 
ground state energy. This happens around x = —0.5 
where Xf ( a measure of kinetic energy) acquires the max- 
imum amplitude. Eventually, when the number of elec- 
trons reaches 2 per site, electrons are frozen at each sites 
and they cannot hop any more and the ground state en- 



ergy is the classical potential energy (3.6£ in our case). 
Therefore, our slave-boson mean-field calculations is able 
to replicate the exact ground state energy at two doping 
limits, and we expect it should describe the intermedi- 
ate doping well. We comment our ground state energy 
bears a similar trend to the one in Kotliar-Ruckenstein 
slave-boson mean-field approach. 47 



IV. CONCLUSIONS 

In this paper we have studied the Kane-Mele-Hubbard 
model on the honeycomb lattice via the slave-boson 
mean-field approach. We have considered both the case 
of half-filling, which has been addressed earlier in the lit- 
erature via other methods, and the case of doping, which 
has not been previously investigated to the best of our 
knowledge. Our main results are summarized in Fig. 1 
and Fig. 8. 

At half-filling, the slave-boson mean-field approach 
fails to capture the correct physics of weak interactions 
and predicts a SC state (rather than a TBI), but we find 
it correctly predicts a spin liquid phase for intermediate 
interactions and small spin-orbit coupling. This is one 
of the least well understood regimes, and in the pres- 
ence of strong spin-orbit coupling is likely to lead to new 
phases. 15,20,31 It would be interesting to consider models 
with further range interactions (first or second-neighbor) 
to see if they might favor any new phases in the phase 
diagram, and possibly other mean-field ansatz for the 
present case as well. 

With finite doping, the spin liquid and dimerized 
phases become superconducting states characterized by 
finite singlet parings (and the superconducting state at 
half- filling remains a superconducting state). Unfortu- 
nately, all the superconducting states we find do not 



possess any obvious form of topological order. Thus, 
our work leaves largely open the question of how likely 
doping-induced topological superconducting states are to 
appear in models that support interacting topological in- 
sulators at half-filling. We hope our work will help to 
stimulate future studies on the effects of doping topolog- 
ical insulators, including those with longer-range inter- 
actions. Doping three dimensional multi-orbital models 
also seems a promising direction. 15,20 
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Appendix A: Slave-boson self-consistency equations 

In this section, we provide some details on the mean- 
field self-consistency equations and order parameters for 
readers interested in the details of our calculations. To 
obtain the self-consistency equations, we start with the 
ground state energy per unit cell E g = Ef+Eb+E c where 
Ef (Eb) is the ground state energy from the fermionic 
(bosonic) part and E c is an energy constant. We have 



E 



l N z^ 



Ai -2y/A 



1/2 



A x + 2V A< 



where A\ and A2 are defined as 



1/2 



(Al) 



A x ee A 2 + \g\^Ai + gfKXi + \ a \H\i + rf Ag X? 



(A2) 



A 2 ee \gf\h\l + \g\h 2 g 2 A' 2 Xl xl - 2\g\ 2 t 2 g 2 A b A' b Xl Xbx'b + W^loX? + |. 9 | 2 t 2 5 2 A 2 A 2 ^ 2 . (A3) 

The bosonic ground state energy is 

JN E [ V (_2A + U}2 - 4|3<A/ " & A 'f X so\ 2 + v (" 2A + U ^ - 4 '^ A / + 3 2 A' / A so 



-U 



(A4) 



where we have chosen excitation spectra 49 that may give rise to Bose-Einstein condensation at k = 0. The energy 

E c = 2A + 2xfi - 6(d B h A + d A h B )tA f + \2{d A h A + d B h B )A / f X so + 6(d A d B - h A h B )t Xf + 6< (A h A/ + XbXf) 

+ K + h 2 B ) (-A + n + GXsox'j) - (4 + d 2 B )(x-U + » + 6A S ox'/) - 12A so (A' 6 A^ + X ' b x'f)- (A5) 

Taking the derivative of E g with respect to the order parameters, we immediately obtain the self-consistency 
equations: 



1 



(d 2 A + d 2 B -h 2 A -h 2 B ), 



(A6) 



Xb = h A h B - d A d B , 



(A7) 



Xb = -^{d 2 A + d 2 B ~h 2 A - h%), 



(A8) 



d A (\ -U + n + GXsox'f) + 3h B tA f - 6h A A' f X so - 3d B tx/ = 0, 

d B (X - U + /i + 6A S oX/) + 3h A tA f - 6h B A' f X so - 3d A t Xf = 0, 

h A (-X + fi + 6AsoX/) - 'id B tA f + 6d A A' f X so - 3h B txt = 0, 

h B (-X + n + 6X so x'f) - 3d A tA f + 6d s A}A so - 3h A t Xf = 0, 



(A9) 
(A10) 
(All) 
(A12) 



XI 



Nt t—' 



6Nt 



2|<?|^X6-g)M 



2,51 * Xb + t ]/d2 



(A13) 



A f = 



-E 

N't ^ 



6Nt 



2| 5 | 2 t 2 A 6 + A 2 |JM 



2\g\VAb-Xl f 5 )/d 2 



(A14) 



A b = d A h B + h A d B 



A' 
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"JnY. \9\ ( (Iff \ tA f + 92 A' f A so ) /d 3 + (l#A/ - 32 A}A so ) /d 4 ) 


(A15) 


Aso v^ 

12N ^ 

V 


;(«*-7r)/* + (*!* + Tr>* 


1 


(A16) 



A' b = h A d A + h B d B - 



6N ^ 

k 



g 2 ( - (\g\tA f +g 2 A' f X so ) /d 3 + (\g\tA f - g 2 A}A S o) M 



(A17) 



Xf 



12iV 



E 



(2^' b -j 8 )/d 1 + (2^ + j 8 )M 

"5 «5 



(A18) 



2 = d A + d B + h A + h B 



l -y 



2N 



i4 '' ^t4 + £) +2 (- 2A + S / '''- 2 ( 2A + t l '' 



(AW) 



where the d; are defined as follows: 



di, 2 = 2 



A 2 + M 2 t 2 (Ag + ^) +ff ?A| (A' 6 2 + xf) 



T 2 J|#t 2 (A 2 + <? 2 A h 2 A 2 ) x 2 ~ 2\g\H^gfA b A' b Xl oXbX ' b + g\ (A 2 + |#t 2 A 2 ) A 2 oX ' b 2 



(A20) 



rf.3,4 = J(-2A + <7) 2 - \\g\H*b? } T Slsl^A/A^Aso - 4 52 2 A^ 2 A 2 , (A21) 



d 5 = y/\g\*1* (A 2 + 5 2 Af A 2 ) x 2 - 2| 5 ! 2 t 2 5?A 6 A' fe Al X6X / 6 + Si (A 2 + M 2 i 2 A 2 ) A§ x b 2 , (A22) 

d 6 = 2\g\H 2 [(A 2 + g\ A 6 2 A 2 ) Xb - . 9 2 A 6 A b A 2 oXb ] , (A23) 



and 



d 7 = 2\g\ 2 t 2 g\x\ (A' bXb - A bX ' b ) , (A24) 



d 8 = 2g\ (-| 5 | 2 i 2 A h A' bXh + A 2 x ' b + |s| 2 i 2 A 2 x ' b ) , (A25) 



d 9 = 2A(| 5 | 2 i 2 x?+ 5 2 Alox' b 2 ), (A26) 



dio = 2| 5 | 2 i 2 <7 2 xt (AUfc ~ A fcXb ) . (A27) 
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